home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / dassl / ddastp.f < prev    next >
Text File  |  1996-07-19  |  18KB  |  613 lines
  1.       SUBROUTINE DDASTP (X, Y, YPRIME, NEQ, RES, JAC, H, WT, JSTART,
  2.      +   IDID, RPAR, IPAR, PHI, DELTA, E, WM, IWM, ALPHA, BETA, GAMMA,
  3.      +   PSI, SIGMA, CJ, CJOLD, HOLD, S, HMIN, UROUND, IPHASE, JCALC,
  4.      +   K, KOLD, NS, NONNEG, NTEMP)
  5. C***BEGIN PROLOGUE  DDASTP
  6. C***SUBSIDIARY
  7. C***PURPOSE  Perform one step of the DDASSL integration.
  8. C***LIBRARY   SLATEC (DASSL)
  9. C***TYPE      DOUBLE PRECISION (SDASTP-S, DDASTP-D)
  10. C***AUTHOR  PETZOLD, LINDA R., (LLNL)
  11. C***DESCRIPTION
  12. C-----------------------------------------------------------------------
  13. C     DDASTP SOLVES A SYSTEM OF DIFFERENTIAL/
  14. C     ALGEBRAIC EQUATIONS OF THE FORM
  15. C     G(X,Y,YPRIME) = 0,  FOR ONE STEP (NORMALLY
  16. C     FROM X TO X+H).
  17. C
  18. C     THE METHODS USED ARE MODIFIED DIVIDED
  19. C     DIFFERENCE,FIXED LEADING COEFFICIENT
  20. C     FORMS OF BACKWARD DIFFERENTIATION
  21. C     FORMULAS. THE CODE ADJUSTS THE STEPSIZE
  22. C     AND ORDER TO CONTROL THE LOCAL ERROR PER
  23. C     STEP.
  24. C
  25. C
  26. C     THE PARAMETERS REPRESENT
  27. C     X  --        INDEPENDENT VARIABLE
  28. C     Y  --        SOLUTION VECTOR AT X
  29. C     YPRIME --    DERIVATIVE OF SOLUTION VECTOR
  30. C                  AFTER SUCCESSFUL STEP
  31. C     NEQ --       NUMBER OF EQUATIONS TO BE INTEGRATED
  32. C     RES --       EXTERNAL USER-SUPPLIED SUBROUTINE
  33. C                  TO EVALUATE THE RESIDUAL.  THE CALL IS
  34. C                  CALL RES(X,Y,YPRIME,DELTA,IRES,RPAR,IPAR)
  35. C                  X,Y,YPRIME ARE INPUT.  DELTA IS OUTPUT.
  36. C                  ON INPUT, IRES=0.  RES SHOULD ALTER IRES ONLY
  37. C                  IF IT ENCOUNTERS AN ILLEGAL VALUE OF Y OR A
  38. C                  STOP CONDITION.  SET IRES=-1 IF AN INPUT VALUE
  39. C                  OF Y IS ILLEGAL, AND DDASTP WILL TRY TO SOLVE
  40. C                  THE PROBLEM WITHOUT GETTING IRES = -1.  IF
  41. C                  IRES=-2, DDASTP RETURNS CONTROL TO THE CALLING
  42. C                  PROGRAM WITH IDID = -11.
  43. C     JAC --       EXTERNAL USER-SUPPLIED ROUTINE TO EVALUATE
  44. C                  THE ITERATION MATRIX (THIS IS OPTIONAL)
  45. C                  THE CALL IS OF THE FORM
  46. C                  CALL JAC(X,Y,YPRIME,PD,CJ,RPAR,IPAR)
  47. C                  PD IS THE MATRIX OF PARTIAL DERIVATIVES,
  48. C                  PD=DG/DY+CJ*DG/DYPRIME
  49. C     H --         APPROPRIATE STEP SIZE FOR NEXT STEP.
  50. C                  NORMALLY DETERMINED BY THE CODE
  51. C     WT --        VECTOR OF WEIGHTS FOR ERROR CRITERION.
  52. C     JSTART --    INTEGER VARIABLE SET 0 FOR
  53. C                  FIRST STEP, 1 OTHERWISE.
  54. C     IDID --      COMPLETION CODE WITH THE FOLLOWING MEANINGS:
  55. C                  IDID= 1 -- THE STEP WAS COMPLETED SUCCESSFULLY
  56. C                  IDID=-6 -- THE ERROR TEST FAILED REPEATEDLY
  57. C                  IDID=-7 -- THE CORRECTOR COULD NOT CONVERGE
  58. C                  IDID=-8 -- THE ITERATION MATRIX IS SINGULAR
  59. C                  IDID=-9 -- THE CORRECTOR COULD NOT CONVERGE.
  60. C                             THERE WERE REPEATED ERROR TEST
  61. C                             FAILURES ON THIS STEP.
  62. C                  IDID=-10-- THE CORRECTOR COULD NOT CONVERGE
  63. C                             BECAUSE IRES WAS EQUAL TO MINUS ONE
  64. C                  IDID=-11-- IRES EQUAL TO -2 WAS ENCOUNTERED,
  65. C                             AND CONTROL IS BEING RETURNED TO
  66. C                             THE CALLING PROGRAM
  67. C     RPAR,IPAR -- REAL AND INTEGER PARAMETER ARRAYS THAT
  68. C                  ARE USED FOR COMMUNICATION BETWEEN THE
  69. C                  CALLING PROGRAM AND EXTERNAL USER ROUTINES
  70. C                  THEY ARE NOT ALTERED BY DDASTP
  71. C     PHI --       ARRAY OF DIVIDED DIFFERENCES USED BY
  72. C                  DDASTP. THE LENGTH IS NEQ*(K+1),WHERE
  73. C                  K IS THE MAXIMUM ORDER
  74. C     DELTA,E --   WORK VECTORS FOR DDASTP OF LENGTH NEQ
  75. C     WM,IWM --    REAL AND INTEGER ARRAYS STORING
  76. C                  MATRIX INFORMATION SUCH AS THE MATRIX
  77. C                  OF PARTIAL DERIVATIVES,PERMUTATION
  78. C                  VECTOR,AND VARIOUS OTHER INFORMATION.
  79. C
  80. C     THE OTHER PARAMETERS ARE INFORMATION
  81. C     WHICH IS NEEDED INTERNALLY BY DDASTP TO
  82. C     CONTINUE FROM STEP TO STEP.
  83. C
  84. C-----------------------------------------------------------------------
  85. C***ROUTINES CALLED  DDAJAC, DDANRM, DDASLV, DDATRP
  86. C***REVISION HISTORY  (YYMMDD)
  87. C   830315  DATE WRITTEN
  88. C   901009  Finished conversion to SLATEC 4.0 format (F.N.Fritsch)
  89. C   901019  Merged changes made by C. Ulrich with SLATEC 4.0 format.
  90. C   901026  Added explicit declarations for all variables and minor
  91. C           cosmetic changes to prologue.  (FNF)
  92. C***END PROLOGUE  DDASTP
  93. C
  94.       INTEGER  NEQ, JSTART, IDID, IPAR(*), IWM(*), IPHASE, JCALC, K,
  95.      *   KOLD, NS, NONNEG, NTEMP
  96.       DOUBLE PRECISION
  97.      *   X, Y(*), YPRIME(*), H, WT(*), RPAR(*), PHI(NEQ,*), DELTA(*),
  98.      *   E(*), WM(*), ALPHA(*), BETA(*), GAMMA(*), PSI(*), SIGMA(*), CJ,
  99.      *   CJOLD, HOLD, S, HMIN, UROUND
  100.       EXTERNAL  RES, JAC
  101. C
  102.       EXTERNAL  DDAJAC, DDANRM, DDASLV, DDATRP
  103.       DOUBLE PRECISION  DDANRM
  104. C
  105.       INTEGER  I, IER, IRES, J, J1, KDIFF, KM1, KNEW, KP1, KP2, LCTF,
  106.      *   LETF, LMXORD, LNJE, LNRE, LNST, M, MAXIT, NCF, NEF, NSF, NSP1
  107.       DOUBLE PRECISION
  108.      *   ALPHA0, ALPHAS, CJLAST, CK, DELNRM, ENORM, ERK, ERKM1,
  109.      *   ERKM2, ERKP1, ERR, EST, HNEW, OLDNRM, PNORM, R, RATE, TEMP1,
  110.      *   TEMP2, TERK, TERKM1, TERKM2, TERKP1, XOLD, XRATE
  111.       LOGICAL  CONVGD
  112. C
  113.       PARAMETER (LMXORD=3)
  114.       PARAMETER (LNST=11)
  115.       PARAMETER (LNRE=12)
  116.       PARAMETER (LNJE=13)
  117.       PARAMETER (LETF=14)
  118.       PARAMETER (LCTF=15)
  119. C
  120.       DATA MAXIT/4/
  121.       DATA XRATE/0.25D0/
  122. C
  123. C
  124. C
  125. C
  126. C
  127. C-----------------------------------------------------------------------
  128. C     BLOCK 1.
  129. C     INITIALIZE. ON THE FIRST CALL,SET
  130. C     THE ORDER TO 1 AND INITIALIZE
  131. C     OTHER VARIABLES.
  132. C-----------------------------------------------------------------------
  133. C
  134. C     INITIALIZATIONS FOR ALL CALLS
  135. C***FIRST EXECUTABLE STATEMENT  DDASTP
  136.       IDID=1
  137.       XOLD=X
  138.       NCF=0
  139.       NSF=0
  140.       NEF=0
  141.       IF(JSTART .NE. 0) GO TO 120
  142. C
  143. C     IF THIS IS THE FIRST STEP,PERFORM
  144. C     OTHER INITIALIZATIONS
  145.       IWM(LETF) = 0
  146.       IWM(LCTF) = 0
  147.       K=1
  148.       KOLD=0
  149.       HOLD=0.0D0
  150.       JSTART=1
  151.       PSI(1)=H
  152.       CJOLD = 1.0D0/H
  153.       CJ = CJOLD
  154.       S = 100.D0
  155.       JCALC = -1
  156.       DELNRM=1.0D0
  157.       IPHASE = 0
  158.       NS=0
  159. 120   CONTINUE
  160. C
  161. C
  162. C
  163. C
  164. C
  165. C-----------------------------------------------------------------------
  166. C     BLOCK 2
  167. C     COMPUTE COEFFICIENTS OF FORMULAS FOR
  168. C     THIS STEP.
  169. C-----------------------------------------------------------------------
  170. 200   CONTINUE
  171.       KP1=K+1
  172.       KP2=K+2
  173.       KM1=K-1
  174.       XOLD=X
  175.       IF(H.NE.HOLD.OR.K .NE. KOLD) NS = 0
  176.       NS=MIN(NS+1,KOLD+2)
  177.       NSP1=NS+1
  178.       IF(KP1 .LT. NS)GO TO 230
  179. C
  180.       BETA(1)=1.0D0
  181.       ALPHA(1)=1.0D0
  182.       TEMP1=H
  183.       GAMMA(1)=0.0D0
  184.       SIGMA(1)=1.0D0
  185.       DO 210 I=2,KP1
  186.          TEMP2=PSI(I-1)
  187.          PSI(I-1)=TEMP1
  188.          BETA(I)=BETA(I-1)*PSI(I-1)/TEMP2
  189.          TEMP1=TEMP2+H
  190.          ALPHA(I)=H/TEMP1
  191.          SIGMA(I)=(I-1)*SIGMA(I-1)*ALPHA(I)
  192.          GAMMA(I)=GAMMA(I-1)+ALPHA(I-1)/H
  193. 210      CONTINUE
  194.       PSI(KP1)=TEMP1
  195. 230   CONTINUE
  196. C
  197. C     COMPUTE ALPHAS, ALPHA0
  198.       ALPHAS = 0.0D0
  199.       ALPHA0 = 0.0D0
  200.       DO 240 I = 1,K
  201.         ALPHAS = ALPHAS - 1.0D0/I
  202.         ALPHA0 = ALPHA0 - ALPHA(I)
  203. 240     CONTINUE
  204. C
  205. C     COMPUTE LEADING COEFFICIENT CJ
  206.       CJLAST = CJ
  207.       CJ = -ALPHAS/H
  208. C
  209. C     COMPUTE VARIABLE STEPSIZE ERROR COEFFICIENT CK
  210.       CK = ABS(ALPHA(KP1) + ALPHAS - ALPHA0)
  211.       CK = MAX(CK,ALPHA(KP1))
  212. C
  213. C     DECIDE WHETHER NEW JACOBIAN IS NEEDED
  214.       TEMP1 = (1.0D0 - XRATE)/(1.0D0 + XRATE)
  215.       TEMP2 = 1.0D0/TEMP1
  216.       IF (CJ/CJOLD .LT. TEMP1 .OR. CJ/CJOLD .GT. TEMP2) JCALC = -1
  217.       IF (CJ .NE. CJLAST) S = 100.D0
  218. C
  219. C     CHANGE PHI TO PHI STAR
  220.       IF(KP1 .LT. NSP1) GO TO 280
  221.       DO 270 J=NSP1,KP1
  222.          DO 260 I=1,NEQ
  223. 260         PHI(I,J)=BETA(J)*PHI(I,J)
  224. 270      CONTINUE
  225. 280   CONTINUE
  226. C
  227. C     UPDATE TIME
  228.       X=X+H
  229. C
  230. C
  231. C
  232. C
  233. C
  234. C-----------------------------------------------------------------------
  235. C     BLOCK 3
  236. C     PREDICT THE SOLUTION AND DERIVATIVE,
  237. C     AND SOLVE THE CORRECTOR EQUATION
  238. C-----------------------------------------------------------------------
  239. C
  240. C     FIRST,PREDICT THE SOLUTION AND DERIVATIVE
  241. 300   CONTINUE
  242.       DO 310 I=1,NEQ
  243.          Y(I)=PHI(I,1)
  244. 310      YPRIME(I)=0.0D0
  245.       DO 330 J=2,KP1
  246.          DO 320 I=1,NEQ
  247.             Y(I)=Y(I)+PHI(I,J)
  248. 320         YPRIME(I)=YPRIME(I)+GAMMA(J)*PHI(I,J)
  249. 330   CONTINUE
  250.       PNORM = DDANRM (NEQ,Y,WT,RPAR,IPAR)
  251. C
  252. C
  253. C
  254. C     SOLVE THE CORRECTOR EQUATION USING A
  255. C     MODIFIED NEWTON SCHEME.
  256.       CONVGD= .TRUE.
  257.       M=0
  258.       IWM(LNRE)=IWM(LNRE)+1
  259.       IRES = 0
  260.       CALL RES(X,Y,YPRIME,DELTA,IRES,RPAR,IPAR)
  261.       IF (IRES .LT. 0) GO TO 380
  262. C
  263. C
  264. C     IF INDICATED,REEVALUATE THE
  265. C     ITERATION MATRIX PD = DG/DY + CJ*DG/DYPRIME
  266. C     (WHERE G(X,Y,YPRIME)=0). SET
  267. C     JCALC TO 0 AS AN INDICATOR THAT
  268. C     THIS HAS BEEN DONE.
  269.       IF(JCALC .NE. -1)GO TO 340
  270.       IWM(LNJE)=IWM(LNJE)+1
  271.       JCALC=0
  272.       CALL DDAJAC(NEQ,X,Y,YPRIME,DELTA,CJ,H,
  273.      * IER,WT,E,WM,IWM,RES,IRES,UROUND,JAC,RPAR,
  274.      * IPAR,NTEMP)
  275.       CJOLD=CJ
  276.       S = 100.D0
  277.       IF (IRES .LT. 0) GO TO 380
  278.       IF(IER .NE. 0)GO TO 380
  279.       NSF=0
  280. C
  281. C
  282. C     INITIALIZE THE ERROR ACCUMULATION VECTOR E.
  283. 340   CONTINUE
  284.       DO 345 I=1,NEQ
  285. 345      E(I)=0.0D0
  286. C
  287. C
  288. C     CORRECTOR LOOP.
  289. 350   CONTINUE
  290. C
  291. C     MULTIPLY RESIDUAL BY TEMP1 TO ACCELERATE CONVERGENCE
  292.       TEMP1 = 2.0D0/(1.0D0 + CJ/CJOLD)
  293.       DO 355 I = 1,NEQ
  294. 355     DELTA(I) = DELTA(I) * TEMP1
  295. C
  296. C     COMPUTE A NEW ITERATE (BACK-SUBSTITUTION).
  297. C     STORE THE CORRECTION IN DELTA.
  298.       CALL DDASLV(NEQ,DELTA,WM,IWM)
  299. C
  300. C     UPDATE Y,E,AND YPRIME
  301.       DO 360 I=1,NEQ
  302.          Y(I)=Y(I)-DELTA(I)
  303.          E(I)=E(I)-DELTA(I)
  304. 360      YPRIME(I)=YPRIME(I)-CJ*DELTA(I)
  305. C
  306. C     TEST FOR CONVERGENCE OF THE ITERATION
  307.       DELNRM=DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
  308.       IF (DELNRM .LE. 100.D0*UROUND*PNORM) GO TO 375
  309.       IF (M .GT. 0) GO TO 365
  310.          OLDNRM = DELNRM
  311.          GO TO 367
  312. 365   RATE = (DELNRM/OLDNRM)**(1.0D0/M)
  313.       IF (RATE .GT. 0.90D0) GO TO 370
  314.       S = RATE/(1.0D0 - RATE)
  315. 367   IF (S*DELNRM .LE. 0.33D0) GO TO 375
  316. C
  317. C     THE CORRECTOR HAS NOT YET CONVERGED.
  318. C     UPDATE M AND TEST WHETHER THE
  319. C     MAXIMUM NUMBER OF ITERATIONS HAVE
  320. C     BEEN TRIED.
  321.       M=M+1
  322.       IF(M.GE.MAXIT)GO TO 370
  323. C
  324. C     EVALUATE THE RESIDUAL
  325. C     AND GO BACK TO DO ANOTHER ITERATION
  326.       IWM(LNRE)=IWM(LNRE)+1
  327.       IRES = 0
  328.       CALL RES(X,Y,YPRIME,DELTA,IRES,
  329.      *  RPAR,IPAR)
  330.       IF (IRES .LT. 0) GO TO 380
  331.       GO TO 350
  332. C
  333. C
  334. C     THE CORRECTOR FAILED TO CONVERGE IN MAXIT
  335. C     ITERATIONS. IF THE ITERATION MATRIX
  336. C     IS NOT CURRENT,RE-DO THE STEP WITH
  337. C     A NEW ITERATION MATRIX.
  338. 370   CONTINUE
  339.       IF(JCALC.EQ.0)GO TO 380
  340.       JCALC=-1
  341.       GO TO 300
  342. C
  343. C
  344. C     THE ITERATION HAS CONVERGED.  IF NONNEGATIVITY OF SOLUTION IS
  345. C     REQUIRED, SET THE SOLUTION NONNEGATIVE, IF THE PERTURBATION
  346. C     TO DO IT IS SMALL ENOUGH.  IF THE CHANGE IS TOO LARGE, THEN
  347. C     CONSIDER THE CORRECTOR ITERATION TO HAVE FAILED.
  348. 375   IF(NONNEG .EQ. 0) GO TO 390
  349.       DO 377 I = 1,NEQ
  350. 377      DELTA(I) = MIN(Y(I),0.0D0)
  351.       DELNRM = DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
  352.       IF(DELNRM .GT. 0.33D0) GO TO 380
  353.       DO 378 I = 1,NEQ
  354. 378      E(I) = E(I) - DELTA(I)
  355.       GO TO 390
  356. C
  357. C
  358. C     EXITS FROM BLOCK 3
  359. C     NO CONVERGENCE WITH CURRENT ITERATION
  360. C     MATRIX,OR SINGULAR ITERATION MATRIX
  361. 380   CONVGD= .FALSE.
  362. 390   JCALC = 1
  363.       IF(.NOT.CONVGD)GO TO 600
  364. C
  365. C
  366. C
  367. C
  368. C
  369. C-----------------------------------------------------------------------
  370. C     BLOCK 4
  371. C     ESTIMATE THE ERRORS AT ORDERS K,K-1,K-2
  372. C     AS IF CONSTANT STEPSIZE WAS USED. ESTIMATE
  373. C     THE LOCAL ERROR AT ORDER K AND TEST
  374. C     WHETHER THE CURRENT STEP IS SUCCESSFUL.
  375. C-----------------------------------------------------------------------
  376. C
  377. C     ESTIMATE ERRORS AT ORDERS K,K-1,K-2
  378.       ENORM = DDANRM(NEQ,E,WT,RPAR,IPAR)
  379.       ERK = SIGMA(K+1)*ENORM
  380.       TERK = (K+1)*ERK
  381.       EST = ERK
  382.       KNEW=K
  383.       IF(K .EQ. 1)GO TO 430
  384.       DO 405 I = 1,NEQ
  385. 405     DELTA(I) = PHI(I,KP1) + E(I)
  386.       ERKM1=SIGMA(K)*DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
  387.       TERKM1 = K*ERKM1
  388.       IF(K .GT. 2)GO TO 410
  389.       IF(TERKM1 .LE. 0.5D0*TERK)GO TO 420
  390.       GO TO 430
  391. 410   CONTINUE
  392.       DO 415 I = 1,NEQ
  393. 415     DELTA(I) = PHI(I,K) + DELTA(I)
  394.       ERKM2=SIGMA(K-1)*DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
  395.       TERKM2 = (K-1)*ERKM2
  396.       IF(MAX(TERKM1,TERKM2).GT.TERK)GO TO 430
  397. C     LOWER THE ORDER
  398. 420   CONTINUE
  399.       KNEW=K-1
  400.       EST = ERKM1
  401. C
  402. C
  403. C     CALCULATE THE LOCAL ERROR FOR THE CURRENT STEP
  404. C     TO SEE IF THE STEP WAS SUCCESSFUL
  405. 430   CONTINUE
  406.       ERR = CK * ENORM
  407.       IF(ERR .GT. 1.0D0)GO TO 600
  408. C
  409. C
  410. C
  411. C
  412. C
  413. C-----------------------------------------------------------------------
  414. C     BLOCK 5
  415. C     THE STEP IS SUCCESSFUL. DETERMINE
  416. C     THE BEST ORDER AND STEPSIZE FOR
  417. C     THE NEXT STEP. UPDATE THE DIFFERENCES
  418. C     FOR THE NEXT STEP.
  419. C-----------------------------------------------------------------------
  420.       IDID=1
  421.       IWM(LNST)=IWM(LNST)+1
  422.       KDIFF=K-KOLD
  423.       KOLD=K
  424.       HOLD=H
  425. C
  426. C
  427. C     ESTIMATE THE ERROR AT ORDER K+1 UNLESS:
  428. C        ALREADY DECIDED TO LOWER ORDER, OR
  429. C        ALREADY USING MAXIMUM ORDER, OR
  430. C        STEPSIZE NOT CONSTANT, OR
  431. C        ORDER RAISED IN PREVIOUS STEP
  432.       IF(KNEW.EQ.KM1.OR.K.EQ.IWM(LMXORD))IPHASE=1
  433.       IF(IPHASE .EQ. 0)GO TO 545
  434.       IF(KNEW.EQ.KM1)GO TO 540
  435.       IF(K.EQ.IWM(LMXORD)) GO TO 550
  436.       IF(KP1.GE.NS.OR.KDIFF.EQ.1)GO TO 550
  437.       DO 510 I=1,NEQ
  438. 510      DELTA(I)=E(I)-PHI(I,KP2)
  439.       ERKP1 = (1.0D0/(K+2))*DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
  440.       TERKP1 = (K+2)*ERKP1
  441.       IF(K.GT.1)GO TO 520
  442.       IF(TERKP1.GE.0.5D0*TERK)GO TO 550
  443.       GO TO 530
  444. 520   IF(TERKM1.LE.MIN(TERK,TERKP1))GO TO 540
  445.       IF(TERKP1.GE.TERK.OR.K.EQ.IWM(LMXORD))GO TO 550
  446. C
  447. C     RAISE ORDER
  448. 530   K=KP1
  449.       EST = ERKP1
  450.       GO TO 550
  451. C
  452. C     LOWER ORDER
  453. 540   K=KM1
  454.       EST = ERKM1
  455.       GO TO 550
  456. C
  457. C     IF IPHASE = 0, INCREASE ORDER BY ONE AND MULTIPLY STEPSIZE BY
  458. C     FACTOR TWO
  459. 545   K = KP1
  460.       HNEW = H*2.0D0
  461.       H = HNEW
  462.       GO TO 575
  463. C
  464. C
  465. C     DETERMINE THE APPROPRIATE STEPSIZE FOR
  466. C     THE NEXT STEP.
  467. 550   HNEW=H
  468.       TEMP2=K+1
  469.       R=(2.0D0*EST+0.0001D0)**(-1.0D0/TEMP2)
  470.       IF(R .LT. 2.0D0) GO TO 555
  471.       HNEW = 2.0D0*H
  472.       GO TO 560
  473. 555   IF(R .GT. 1.0D0) GO TO 560
  474.       R = MAX(0.5D0,MIN(0.9D0,R))
  475.       HNEW = H*R
  476. 560   H=HNEW
  477. C
  478. C
  479. C     UPDATE DIFFERENCES FOR NEXT STEP
  480. 575   CONTINUE
  481.       IF(KOLD.EQ.IWM(LMXORD))GO TO 585
  482.       DO 580 I=1,NEQ
  483. 580      PHI(I,KP2)=E(I)
  484. 585   CONTINUE
  485.       DO 590 I=1,NEQ
  486. 590      PHI(I,KP1)=PHI(I,KP1)+E(I)
  487.       DO 595 J1=2,KP1
  488.          J=KP1-J1+1
  489.          DO 595 I=1,NEQ
  490. 595      PHI(I,J)=PHI(I,J)+PHI(I,J+1)
  491.       RETURN
  492. C
  493. C
  494. C
  495. C
  496. C
  497. C-----------------------------------------------------------------------
  498. C     BLOCK 6
  499. C     THE STEP IS UNSUCCESSFUL. RESTORE X,PSI,PHI
  500. C     DETERMINE APPROPRIATE STEPSIZE FOR
  501. C     CONTINUING THE INTEGRATION, OR EXIT WITH
  502. C     AN ERROR FLAG IF THERE HAVE BEEN MANY
  503. C     FAILURES.
  504. C-----------------------------------------------------------------------
  505. 600   IPHASE = 1
  506. C
  507. C     RESTORE X,PHI,PSI
  508.       X=XOLD
  509.       IF(KP1.LT.NSP1)GO TO 630
  510.       DO 620 J=NSP1,KP1
  511.          TEMP1=1.0D0/BETA(J)
  512.          DO 610 I=1,NEQ
  513. 610         PHI(I,J)=TEMP1*PHI(I,J)
  514. 620      CONTINUE
  515. 630   CONTINUE
  516.       DO 640 I=2,KP1
  517. 640      PSI(I-1)=PSI(I)-H
  518. C
  519. C
  520. C     TEST WHETHER FAILURE IS DUE TO CORRECTOR ITERATION
  521. C     OR ERROR TEST
  522.       IF(CONVGD)GO TO 660
  523.       IWM(LCTF)=IWM(LCTF)+1
  524. C
  525. C
  526. C     THE NEWTON ITERATION FAILED TO CONVERGE WITH
  527. C     A CURRENT ITERATION MATRIX.  DETERMINE THE CAUSE
  528. C     OF THE FAILURE AND TAKE APPROPRIATE ACTION.
  529.       IF(IER.EQ.0)GO TO 650
  530. C
  531. C     THE ITERATION MATRIX IS SINGULAR. REDUCE
  532. C     THE STEPSIZE BY A FACTOR OF 4. IF
  533. C     THIS HAPPENS THREE TIMES IN A ROW ON
  534. C     THE SAME STEP, RETURN WITH AN ERROR FLAG
  535.       NSF=NSF+1
  536.       R = 0.25D0
  537.       H=H*R
  538.       IF (NSF .LT. 3 .AND. ABS(H) .GE. HMIN) GO TO 690
  539.       IDID=-8
  540.       GO TO 675
  541. C
  542. C
  543. C     THE NEWTON ITERATION FAILED TO CONVERGE FOR A REASON
  544. C     OTHER THAN A SINGULAR ITERATION MATRIX.  IF IRES = -2, THEN
  545. C     RETURN.  OTHERWISE, REDUCE THE STEPSIZE AND TRY AGAIN, UNLESS
  546. C     TOO MANY FAILURES HAVE OCCURED.
  547. 650   CONTINUE
  548.       IF (IRES .GT. -2) GO TO 655
  549.       IDID = -11
  550.       GO TO 675
  551. 655   NCF = NCF + 1
  552.       R = 0.25D0
  553.       H = H*R
  554.       IF (NCF .LT. 10 .AND. ABS(H) .GE. HMIN) GO TO 690
  555.       IDID = -7
  556.       IF (IRES .LT. 0) IDID = -10
  557.       IF (NEF .GE. 3) IDID = -9
  558.       GO TO 675
  559. C
  560. C
  561. C     THE NEWTON SCHEME CONVERGED,AND THE CAUSE
  562. C     OF THE FAILURE WAS THE ERROR ESTIMATE
  563. C     EXCEEDING THE TOLERANCE.
  564. 660   NEF=NEF+1
  565.       IWM(LETF)=IWM(LETF)+1
  566.       IF (NEF .GT. 1) GO TO 665
  567. C
  568. C     ON FIRST ERROR TEST FAILURE, KEEP CURRENT ORDER OR LOWER
  569. C     ORDER BY ONE.  COMPUTE NEW STEPSIZE BASED ON DIFFERENCES
  570. C     OF THE SOLUTION.
  571.       K = KNEW
  572.       TEMP2 = K + 1
  573.       R = 0.90D0*(2.0D0*EST+0.0001D0)**(-1.0D0/TEMP2)
  574.       R = MAX(0.25D0,MIN(0.9D0,R))
  575.       H = H*R
  576.       IF (ABS(H) .GE. HMIN) GO TO 690
  577.       IDID = -6
  578.       GO TO 675
  579. C
  580. C     ON SECOND ERROR TEST FAILURE, USE THE CURRENT ORDER OR
  581. C     DECREASE ORDER BY ONE.  REDUCE THE STEPSIZE BY A FACTOR OF
  582. C     FOUR.
  583. 665   IF (NEF .GT. 2) GO TO 670
  584.       K = KNEW
  585.       H = 0.25D0*H
  586.       IF (ABS(H) .GE. HMIN) GO TO 690
  587.       IDID = -6
  588.       GO TO 675
  589. C
  590. C     ON THIRD AND SUBSEQUENT ERROR TEST FAILURES, SET THE ORDER TO
  591. C     ONE AND REDUCE THE STEPSIZE BY A FACTOR OF FOUR.
  592. 670   K = 1
  593.       H = 0.25D0*H
  594.       IF (ABS(H) .GE. HMIN) GO TO 690
  595.       IDID = -6
  596.       GO TO 675
  597. C
  598. C
  599. C
  600. C
  601. C     FOR ALL CRASHES, RESTORE Y TO ITS LAST VALUE,
  602. C     INTERPOLATE TO FIND YPRIME AT LAST X, AND RETURN
  603. 675   CONTINUE
  604.       CALL DDATRP(X,X,Y,YPRIME,NEQ,K,PHI,PSI)
  605.       RETURN
  606. C
  607. C
  608. C     GO BACK AND TRY THIS STEP AGAIN
  609. 690   GO TO 200
  610. C
  611. C------END OF SUBROUTINE DDASTP------
  612.       END
  613.